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Preface 


In  the  fall  of  1979,  a  program  of  balloon-borne  mass  spectrometric  measure¬ 
ments  was  commenced  at  this  laboratory.  This  work  has  been  continued  up  to  the 
present  time,  hut  little  support  in  the  way  of  theory  has  been  undertaken.  The 
present  work  is  addressed  towards  making  a  start  in  this  direction;  it  deals  with 
the  problems  encountered  in  relating  measured  quantities  to  number  densities  of 
stratospheric  ionic  and  neutral  species. 

I  thank  Pat  Hench  for  programming  the  numerical  solutions  and  Tony  Quesada 
for  bringing  several  useful  references  to  my  attention. 
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I 

Collection  of  Ions  in  the  Stratosphere 


I.  INTROIM  'CTION 

I .  I  Hierarchy  of  Time  ConxIanU 

The  problem  of  sampling  in  the  atmosphere  is  that  of  finding  how  the  number 
densities  and  fluxes  of  species  we  collect  relate  to  the  undisturbed  ambient  quan¬ 
tities.  The  nature  of  the  problem  will  vary  widely  depending  on  the  altitude  re¬ 
gime  considered.  Here  we  consider  collection  using  a  balloon  to  carry  the  exper¬ 
iment  in  the  regime  between  20  and  40  km,  so  that  continuum  conditions  may  be 
taken  to  prevail.  We  envision  the  sampling  to  occur  through  a  small  orifice  but  do 
not  follow  the  flow  through  the  orifice.  Instead,  we  ask  the  question:  In  the  ab¬ 
sence  of  the  orifice,  what  is  the  flux  of  charged  particles,  here  predominantly  pos¬ 
itive  and  negative  ions,  to  the  collecting  surface  in  terms  of  the  ambient  number 
densities  ° 

What  happens  to  particles  on  entry  info  the  detecting  orifice  and  further  into  a 
detecting  instrument  is  a  separate  problem  of  interest  in  its  own  right. 

In  this  report,  we  do  not  undertake  detailed  calculations  of  the  complicated 
geometries  actually  utilized.  Rather,  we  discuss  the  collection  process  in  general 
terms  to  assess  the  relative  importance  of  the  various  processes  involved,  and, 
with  an  eye  to  these  processes,  we  solve  some  relevant  prototypical  problems  for 
simple  geometries. 

(Received  for  publication  9  December  19851 
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There  are,  to  begin  with,  three  transport  processes  involved  in  the  collection 
of  charged  particles  in  the  stratosphere:  bulk  fluid  flow,  diffusion,  and,  if  electric 
fields  are  present,  mobility  effects.  To  estimate  the  relative  importance  of  these 
effects  as  well  as  to  compare  them  with  other  non-transport  processes,  we  de¬ 
velop  a  hierarchy  of  time  constants,  T..  First,  we  deal  with  the  time  constant  for 
bulk  fluid  flow.  This  is  given  by  a  characteristic  length  L,  divided  by  a  veloc¬ 
ity  u  .  Now,  the  balloon  itself  moves  with  the  ambient  wind  and  velocities  relative 
to  the  sampling  platform  are  generated  either  by  horizontal  gradients  or  by  local 
fluctuations  in  wind  speed.  If  we  take  3  <  L  <  100cm  and  10  <  uq  <  100 cm/sec, 
we  have  0.03  <  <  io  sec  where  the  subscript  B  refers  to  bulk  flow. 

For  diffusion,  typical  time  constants  are  ordinarily  taken  to  be 

T  L2  _  3  L2  ' 

D  4D  4  A  v 

where  I)  =  the  diffusion  coefficient,  A  =  the  mean  free  path,  and  v  *  root  mean 
square  thermal  velocity.  In  the  present  instance,  however,  this  formula  may  not 
be  valid.  Bulk  flow  across  a  surface  creates  a  boundary  layer,  and,  in  the  ab¬ 
sence  of  diffusion,  the  number  densities  at  the  free  edge  of  the  boundary  are  the 
undisturbed  values.  The  characteristic  dimension  for  diffusion  for  transport 
across  this  layer  then  becomes  the  boundary  layer  thickness  6  *  L/  yile  where  Re 
is  the  hulk  flow  Reynolds  number.  We  have  then. 


TD  — 


L  2(i 


4  A  v  uc  p  L 


4  A  vRe 

(p  is  the  coefficient  of  viscosity,  p  the  density) 


and  putting  p  *  y  Apv, 


T  =  -1  — 
U  4  u0 


We  thus  see  that  the  range  of  is  of  the  same  order  as  that  of  T^.  This  is  an 
interesting  result.  It  implies  that,  although  isolated  diffusive  times  may  be  long 
compared  to  isolated  bulk  flow  times,  once  bulk  flow  occurs,  it  enhances  diffusion 
speeds  so  that  these  latter  "keep  step"  with  flow,  as  far  as  order  of  magnitude  of 
influence  is  concerned;  and  further,  that  this  occurs  independently  of  altitude.  If 
this  interpretation  is  correct,  great  caution  must  be  used  before  neglecting  diffu¬ 
sive  effects  as  being  "small"  in  the  sampling  process. 

We  next  consider  the  time  constants  associated  with  mobility  phenomena.  The 


ratio  of  mobility  times  to  diffusion  times  is  given  by 


Vfe)  (?)“(£)(*?) 


Av 


Re  =  mv  Re  *  rrl^  Re 

eX0o  e0o 


Here  K.  E,  0o>  e,  and  m  are  the  mobility  coefficient,  electric  field,  electric  po¬ 
tential.  electronic  charge,  and  particle  mass,  respectively.  We  take  0  to  lie 
between  1  and  20  V  (typical  collecting  voltages),  and  take  L  =  10cm  and  v  = 

30  cm /sec.  We  have 


Re  = 


Pv  oL 


p  300 

Pn  o7l5 


P 

2000  — 


where 


is  ihe  density  at  sea  level,  and 


Alt  (km) 

PlPo 

Re 

mv  /e  (V) 

0(V) 

tM/tD 

20 

1/20 

100 

1/25 

1 

4.0 

40 

1/400 

5 

1/25 

20 

0.01 

Thus,  the  range  of  overlaps  the  range  of  Tg.  This  completes  the  estimates  of 
transport  times. 

Finally,  we  turn  to  the  question  of  generation  and  recombination  times.  Since 
the  charged  particle  number  density  no>  the  generation  rate  q,  and  the  recombina¬ 
tion  rate  a  are  connected  by  n„  =  Vq  /a  ,  the  generation  and  recombination  times 

v  -1/2 

are  identical  and  are  given  by  tg  =  Tp  =  <<*  q )  .  We  then  have  the  following: 


Alt  (km) 

,  -3  -1. 

q  (cm  sec  ) 

,  3  -1. 

a  (cm  sec  ) 

nQ(cm  3) 

TG'TR 

20 

15 

5  x  10'7 

~5500 

360 

40 

0.  75 

5  x  10-8 

~  4000 

5200 

Thus,  T, 


G’  tr 


>  >  rB.  rD.  rM 


We  are  now  in  position  to  state  that  when  we  collect  ions  in  the  stratosphere, 
the  effects  of  bulk  flow,  diffusion,  and  mobility  tend  to  dominate  those  of  genera¬ 
tion  and  recombination.  Although  the  latter  two  establish  the  equilibrium  number 
densities,  it  will  often  be  possible  to  fix  nQ  as  a  boundary  condition  in  the  collec¬ 
tion  process. 


1.2  Space  Charge 

In  considering  the  effects  of  space  charge,  it  is  expedient  to  think  in  terms 
other  than  those  of  time  constants. 

Any  biased  collecting  system  (probe)  placed  in  a  plasma  generates  space 
charge  (that  is.  separation  of  charges).  If,  however,  the  distribution  of  charge  is 
such  that  the  changes  induced  in  the  primary  field  are  small,  the  primary  field 
may  be  used  to  obtain  collected  currents,  a  much  simpler  problem  than  the  self- 
consistent  solution  for  fields  and  fluxes.  The  criterion  as  to  when  such  approxi¬ 
mation  is  permissible  is,  in  general,  connected  with  the  Debye  length  kD.  If 
Xp  <  L,  then  space  charge  effects  tend  to  be  either  dominant  or  at  least  substan¬ 
tial;  while  if  Ad>>  L,  they  are  small  or  negligible.  In  the  present  instance,  the 
Debye  length  for  the  stratosphere  is  of  the  order  of  2  cm,  while  material  dimen¬ 
sions,  L,  are  such  that  10  <  L  <  100  cm.  The  indications  are  thus  that  space 
j  charge  effects  will  be  substantial. 

)  It  is  known,  however,  that  the  ratio  of  Debye  length  to  material  dimension  is 

not  the  only  dimensionless  parameter  related  to  space  charge  effects.  A  second 
and  independent  characteristic  length  Ag  *  e.  o/nQLe,  where  eq  is  the  permitti¬ 
vity  of  free  space,  may  also  be  formed.  The  independence  of  Ag  from  AD  is  obvi¬ 
ous:  Aq  contains  only  parameters  related  to  plasma  properties,  while  both  and 
L  are  boundary  properties. 

The  ratio  yg/L  is  the  reciprocal  of  the  ratio  of  secondary  (that  is,  space 
charge  generated)  electric  fields  to  primary  electric  field.  Thus,  for  the  second¬ 
ary  field.  Eg,  we  have  on  dimensional  grounds. 


E 


S 


e  o 


and 


2 

Eg  n0eL  L 
Eo  Eo*o  LS 


This  dimensionless  ratio  may  be  expressed  in  terms  commonly  appearing  in  probe 
theory 

v4- 

E0  A2d  V  kT/ 


where  k  is  Boltzmann's  constant,  and  T  is  the  temperature. 
For  nQ  »  5  x  109/m3, 

4 


Taking  L  >  lm  and  *  IV, 


5a 

E_ 


=  80 


■'MAX 


while  for  L  *  0. 1  m  and  $  *  20  V 

o 


E, 


=  0.04 


'MIN 


This  seems  to  indicate  that  under  some  conditions  space  charge  will  not  play 
an  important  role  in  collection.  This  conclusion  is  the  opposite  of  that  drawn  on 
the  basis  of  Debye  length  considerations. 

Although  resolution  of  the  role  of  space  charge  is  not  a  minor  matter,  includ¬ 
ing  space  charge  would  lead  into  the  domain  of  continuum  probe  theory,  an  en¬ 
deavor  that  cannot  be  pursued  here.  We  confine  our  calculations  to  those  that  ig¬ 
nore  space  charge. 

1.3  (iriimtin 

The  overall  geometry  of  the  balloon  platform  and  of  the  environs  of  the  mass 
spectrometer  entrance  aperture  are  so  complex  that  a  very  detailed  numerical 
calculation  would  be  required  to  treat  the  full  surface  boundary  properly.  Here, 
we  consider  two  highly  simplified  aspects  of  the  collection  problem.  In  Section  3, 
we  consider  a  possible  effect  on  collection  of  the  presence  of  the  whole  platform 
taken  to  be  a  sphere,  in  the  presence  of  diffusion,  generation,  and  recombination 
of  charged  particles.  In  Section  2,  we  consider  a  two-dimensional  rectangular 
geometry  in  the  immediate  vicinity  of  the  instrument  entrance  aperture.  There  is 
a  possible  application  for  this  model  calculation.  This  is  an  attempt  to  gain  con¬ 
trol  of  the  sampling  by  creating  a  bulk  flow  sufficient  to  dominate  random  fluctu¬ 
ations  present  in  ambient  winds  (Figure  1-1).  The  idea  is  that,  bv  controlling  the 
location  of  the  housing  and  the  speed  of  the  bulk  flow  induced  by  the  fan  as  well  as 
the  electric  field  intensity,  the  collection  rate  may  be  made  independent  of  the 
ambient  winds.  At  the  same  time,  again  by  proper  control  of  entry  conditions,  it 
should  be  possible  to  avoid  contamination  arising  in  the  control  surface  itself. 
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Figure  1-1.  Schematic  Representation  of  Simplified  Collecting  System 


The  calculations  done  in  Section  2  are  preliminary  to  a  model  for  such  a  collection 
system. 

2.  MFFI'SIOV  MOBILITY.  \\l>  FLOW  WITH  RECTANGl  I.AR  GEOMETRY 
2.1  Series  Solution 

In  this  section,  we  solve  a  simple  boundary  value  problem  involving  diffusion, 
bulk  flow,  and  mobility.  Under  conditions  for  which  a  definite  known  mobility  law 
is  valid  and  there  is  also  no  space  charge,  there  is  then  no  real  distinction  be¬ 
tween  bulk  and  mobility  flows.  The  electric  field  is  solved  for  as  a  separate  prob¬ 
lem,  and  the  total  flow  field  resulting  from  both  electrical  and  mechanical  forces 
is  then  a  known  function  that  is  inserted  into  the  continuity  equation. 

One  of  the  purposes  of  this  calculation  is  to  illustrate  the  role  that  diffusion 
may  play  in  modifying  both  density  distributions  and  fluxes.  It  is  often  taken  for 
granted  that,  in  the  presence  of  substantial  flow  fields,  diffusion  effects  may  be 
neglected.  This  is  universally  done,  for  example,  in  the  discussion  of  Gerdien 


condensers.  In  fact,  diffusion  may  play  a  significant  role  in  ion  collection  even 
when  the  pertinent  dimensionless  parameters  indicate  that  flow  dominates  diffusion. 

A  further  reason  for  including  diffusion  in  sampling  calculations  is  to  obtain 
consistency  in  the  application  of  boundary  conditions  to  the  relevant  differential 
equations.  Diffusion  equations  are  of  higher  order  than  inviscid  flow  equations 
and  thus  permit  the  specification  of  ion  number  densities  on  both  collecting  and 
other  surfaces  (or  at  infinity).  This  cannot  be  done  when  diffusion  is  neglected. 

The  boundary  value  problem  we  investigate  here  will  be  solved  in  three  ways: 
(1)  by  a  series  solution;  (2)  by  a  numerical  machine  solution;  and  (3)  by  a  singular 
perturbation  method.  These  three  methods  are  complementary  and  contribute  to 
the  evaluation  of  validity  of  numerical  solutions,  a  benefit  in  extending  the  solu¬ 
tions  to  situations  that  cannot  be  solved  analytically. 

The  continuity  equation  for  number  density  of  ions,  n',  undergoing  diffusion, 
bulk  flow,  and  mobility  effects  is 

V- [n'  <u  +  fi~E)  -  DVn']  =  0  (2-1) 

where  tT  is  the  bulk  flow  velocity,  and,  in  this  section,  fl  is  the  mobility  coeffi¬ 
cient,  E  the  electric  field,  and  D  the  diffusion  coefficient.  If  we  assume  that  D  is 
constant,  we  have 


n ' V •  (u  +  fi E>  +  (u  +  fjE)  •  */n'  -  DV~  n'  =  0  (2-2) 

For  incompressible  flow  and  no  space  charge,  V*  u  =  V*  E  =  0.  Further,  we  let 

tT  =  u  (y)  and  fjE  -  ~v  -  v  (see  Figure  1-1)  so  that 
*  y 


DV2n'  -  u  ££  -  v  dn' 


dx' 


dy< 


=  0 


(2-3) 


Now,  letting  the  separation  between  the  housing  plates  be  L,  and  using  rectangular 
coordinates,  set 


y  -  y'/L ,  x  -  x'/L,  n  =  n'/no 


Then 


where 


and  in  this  section,  we  use  a  different  definition  of  p. 


and  n(>  is  the  ambient  ion  density.  The  boundary  conditions  for  this  problem  are: 
n  (x,  0)  *  n  (x,  1)  =  0  (2-5 


n  (0,  y)  =  1 


n  (“•,  y)  =  0 


Now,  let  u  *  u  and  v  *  v  be  constants.  Then  the  Eq.  (2-4)  is  separable,  and, 

tt  y 

for  the  given  boundary  conditions,  the  solution  (details  of  which  are  omitted)  is 
n  =  87Texp/-° - A  £  «?XP  f- (0  2  +  P2  ’  4ff2  j2)  -1  —5-^ - y~2 

\  2  /  j=l  L  2  J  p2  +  Air2  jZ 

*[l  +  (-l)-*  +  1  e"*^2]  sin  »rjy  (2-8) 

To  find  the  flux  to  the  collecting  surface  y  =  0,  we  reintroduce  the  primes  which 
designate  the  original  dimensional  variables.  We  have 


Flux-F'  =  -DV'n’  I  =  -DaTT 


y'  =  O 


dy'  I  , 

*•  y 1  *  o 


D  dn 

I 


Let  F  =  F'/n  v.  Then 
o 


p  -  .  EL  dn  .  1  dn 

LV  dy  P  dy 


(2-10) 


F  is  the  non-dimensional  ratio  of  the  observed  flux  to  the  diffusion-free  flux.  We 


(2-11) 


-  pF  =  8n2  e  0/2  x 


53  exp 


r2  +  p2+4,2j2)1/2f| 


P2  +  4»r2  j2 


*  [l  +  (-l)j  +  1  e'P/2] 


We  note  that,  with  our  choice  of  coordinate  axes  and  flow  directions,  u  and  v 
are  positive  and  negative  respectively;  hence,  in  Eqs.  (2-8)  and  (2-11)  (  and  in 
general  herein),  o  and  p  are  also  positive  and  negative  respectively. 

The  utility  of  Eqs.  (2-8)  and  (2-11)  depends  to  a  large  extent  on  rates  of  con¬ 
vergence,  and  this  latter  clearly  depends,  in  turn,  on  values  of  o  and  p.  If  we 
assume  that  O  fa]  =  O  [-P].  it  is  seen  that  both  the  summed  exponential  and 
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(P  +  4rr'j  )  terms  start  to  decrease  when  Ofj]  =  P/27T.  Thus,  to  obtain  conver¬ 


gence,  we  must  have  OfjMAX]  >  >  pl2n,  where  jMAX  is  the  terminal  value  of 
j  in  the  sum. 

This,  however,  is  only  a  necessary,  not  a  sufficient,  condition  for  these  equa¬ 
tions  to  represent  valid  solutions;  for  we  must  guarantee  not  only  that  the  series 
converges,  but  that,  in  a  practical  sense,  it  converges  to  correct  values.  If  the 
series  is  summed  by  digital  computer,  we  may  expect  approximately  12-place  ac¬ 
curacy.  When  O[o]  =  Of-p]  a  100,  initial  terms  in  the  series  (Ofj]  a  1)  are  of 
order 


e*p{§[l  (1  -  yiljj 


which,  depending  on  the  value  of  x,  may  be  as  great  as  ~10A  .  Since,  for  n,  the 
correct  converged  values  must  be  =  1,  it  is  clear  that  for  such  large  values  of  a 
and  -p,  the  series  sum  will  be  meaningless,  even  though  the  series  converges. 
Actual  summing  bears  these  ideas  out.  At  0  =  -  p  =  80,  the  results  are  without 
meaning;  at  a  =  -p  =  60,  the  results  appear  to  be  valid.  Thus,  a  practical  limit 
is  set  on  the  range  of  values  of  a  and  -p  for  which  this  series  solution  can  be  uti- 


Because  of  the  limitations  on  use  of  the  series  solution,  it  is  important  to 
know  what  range  of  values  of  a  and  p  to  expect  in  practice.  For  this  purpose,  we 
choose  u,  -v  »  10  ml  sec  and  L  =  1  m  as  upper  limits,  and  u,  -v  s  1  m/sec,  1.  - 
0.  1  m  as  lower,  we  thus  obtain  (using  values  of  X  from  Appendix  C) 
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The  range  of  o  thus  includes  both  values  which  can  and  cannot  be  summed  in  se¬ 
ries  form.  This  points  up  the  need  for  alternative  solutions. 

We  now  present  graphically  some  of  the  results  obtained  from  Eqs.  (2-8)  and 
(2-11).  Figure  2-1  shows  fluxes,  -p F,  both  with  and  without  diffusion,  as  a  func¬ 
tion  of  x,  with  a  =  -p  as  a  parameter.  Even  at  a  =  60,  there  are  still  substantial 


x 


Figure  2-1.  Normalized  Flux  -p F  vs  x.  With  -p  *  a  as  a 
Parameter.  Solid  lines  are  obtained  from  series  solution; 
dashed,  solution  without  diffusion 


Influences  due  to  the  presence  of  diffusion.  (In  looking  at  these  graphs,  the  values 
of  flux  near  x  •  0,  showing  an  actual  or  incipient  increase,  need  cause  no  concern. 
This  increase  is  an  artifact  due  to  the  singularity  in  boundary  conditions  at  x  =  y 
*  0. )  Figure  2-2  shows  fluxes  -pF  as  a  function  of  o  for  P  -  -30  and  three  loca¬ 
tions  x  again  with  the  corresponding  diffusionless  cases,  also  shown  in  dashed 
lines.  The  approach  to  the  diffusionless  case  (-pF  *  30  for  x  <  -  o/p  ;  -pF  *  0 
for  x  >  -  o/p)  with  increasing  values  of  o  is  seen. 

Figures  2-3  and  2-4  show  number  densities  as  functions  of  x  and  y,  respec¬ 
tively,  for  o  *  -p  ■  10.0  obtained  by  summing  Eq.  (2-8).  The  corresponding  dif¬ 
fusionless  series  have  not  been  drawn  in,  but,  as  might  be  expected  with  o  =  -p 
=  10,  the  solutions  with  diffusion  differ  very  substantially  from  those  without. 


Figure  2-2.  Normalized  Flux  -  PF  vs  o  .  With  x  as  a  Parameter  and  -p 
=  30.  Solid  lines  are  obtained  from  series  solution;  dashes,  without 
diffusion 
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2.2  Nvnerkal  Solution 


In  this  section,  we  take  up  solutions  of  Eq.  (2-4)  by  numerical  methods.  We 
do  this  for  two  reasons:  first,  as  noted,  to  extend  solutions  to  larger  values  of  a 
and  -p,  not  attainable  by  the  series  solution,  and  second,  to  explore  methods  of 
solution  for  more  complicated  equations  for  which  analytic  solutions  are  not  avail¬ 
able.  Thus,  for  example,  uQ  could  be  a  boundary  layer  flow  and  v  a  velocity  in¬ 
duced  by  an  electrode  with  a  finite  edge  leading  to  non-constant  electric  fields. 
Details  of  the  method  of  solution  are  presented  in  Appendix  A.  As  seen  from  this 
appendix,  the  transformation  from  differential  to  a  difference  equation  is  accom¬ 
plished  in  a  simple,  obvious  way.  Subsequent  tests  made  for  values  of  o  ,  -p  <  60 
indicate  that  the  solution  was  stable  and  converging  to  the  correct  (series)  value; 
however,  the  accuracy  attained  leaves  much  to  be  desired.  Figure  2-5  shows  val- 


NUMBER  OF  GRID  POINTS  k ,1 


Figure  2-5.  Number  Density  n,  vs  Number  of  Grid  Points  k  »  t  for  Four  Selected 
Values  of  (x,  y)  and  o  ■  -p  ■  10.  0  as  Obtained  From  Numerical  Solution 


ues  of  n  obtained  from  the  numerical  solution  for  four  selected  values  (x,  y),  as  a 
function  of  the  number  of  increments  on  a  side  of  the  grid  mesh.  An  equal  number 
of  grid  increments  in  the  x  and  z  directions  are  used.  The  horizontal  lines  repre¬ 
sent  the  series  solution.  The  numerical  values  are  seen  to  be  converging  towards 
the  series  solutions,  but  even  with  k  •  t  ■  10,  the  values  for  x,  y  ■  (0.  21,  0.50), 
(0.  21,  0. 10),  and  (1.02,  0. 10)  differ  from  the  series  values  by  «*  3  percent,  2  per¬ 
cent,  and  6  percent,  respectively.  For  large  values  of  x,  the  agreement  becomes 
rapidly  worse.  Use  of  more  grid  points  is  not  feasible  because  even  the  100-by- 
100  mesh  requires  ~25  minutes  of  central  processor  time  on  our  Cyber  750,  and 
CP  time  increases  as  the  fourth  power  of  the  number  of  grid  increments  on  a  side. 

Some  results  of  a  numerical  solution  for  o  *  -p  *  100.0,  values  not  amenable 
to  solution  in  series  form,  are  shown  in  Figure  2-6.  Here,  n  is  plotted  versus  x, 
with  y  -  0. 5  fixed  and  the  number  of  grid  points  per  side  as  a  parameter.  Also 
plotted  on  the  grid  are  the  results  of  a  perturbation  calculation  (discussed  in  Sec- 


Figure  2-6.  Number  Density  n,  vs  x  for  y  *  0.  5,  o  -  -p  *  100,  With  k  =  1  as  a 
Parameter,  as  Obtained  From  Numerical  Solution.  Dashed  curve  is  the  perturba¬ 
tion  solution 


tion  2.3).  We  do  see,  however,  that  the  numerical  solution  continues  to  approach 
the  perturbation  solution  as  the  number  of  grid  points  is  increased,  an  indication 
of  stability  if  one  assumes  the  perturbation  solution  to  be  valid. 

The  slowness  of  convergence  as  well  as  the  lack  of  accuracy,  however,  are 
not  reassuring.  Before  numerical  methods  such  as  these  can  be  used  to  solve 
flow/diffusion  problems,  it  will  be  necessary  to  investigate  questions  of  accuracy 
and  rapidity  of  convergence  with  more  care. 

2.3  Singular  (Boundary  Layer)  Perturbation  Solution 

For  larger  values  of  o,  p.  we  expect  that  some  sort  of  perturbation  method 

will  be  applicable  in  finding  solutions  of  Eq.  (2-4).  Since  the  large  flow  terms  are 

of  lower  order  derivatives  than  the  smaller  (diffusion)  terms,  this  will  be  a  singu- 

1  2 

lar  perturbation  problem.  '  These  problems  are  not  always  straightforward  to 
solve,  and,  in  the  present  instance,  we  intend  to  make  generous  use  of  intuitive 
arguments  rather  than  rigor  in  obtaining  a  solution. 

Zauderer  (p.  595)  gives  a  general  singular  perturbation  solution  to  the  follow¬ 
ing  problem,  in  which  e  is  a  small  positive  quantity. 

e  (u  +  u  )  +  u  +  bu  =0  in 
xx  yy  x  y 

0  <  x  <  “  and  0  <  y  <  L 

and  b  *  constant,  subject  to 

u(x,  o)  =  f(x);  u(o,  y)  *  g (y ) :  u  <x,  L)  *  h(x). 

The  solution  is 

u(x,  y)  *  h  ^x  -  }  *  £f(x>  -  h^x  +  ^je’(b^e  ,v 

-  g(v)  -  h 

It  is  seen  immediately  that  this  problem  has  much  in  common  with  the  one  we  are 
trying  to  solve.  Both  the  basic  partial  differential  equation,  and,  if  we  let  L  1  1, 

1.  Carrier,  George  F. ,  and  Pearson,  Carl  E.  (1976)  Partial  Differential  Equa¬ 

tions,  Academic  Press,  New  York. 

2.  Zauderer,  Erich,  (1983)  Partial  Differential  Equations  of  Applied  Mathema¬ 

tics.  John  Wiley  &  Sons,  New  York. 


the  region  of  investigation,  are  the  same.  However,  the  boundary  values  in  our 
case,  although  simple,  have  certain  peculiarities  which  make  it  preferable  to  solve 
our  problem  directly  without  recourse  to  the  general  formula.  We  now  proceed  to 
do  this. 

First,  replacing  u  by  n  (in  accord  with  our  previous  notation),  the  bulk  solu¬ 
tion  (that  is,  outside  of  any  boundary  layer)  may  be  written  down  by  inspection 


n  =  J  d(x')  dx' 
y  -  1  +  bx 


(2-12) 


where  b  =  - p/a  >  0  and  fl(x')  is  the  Dirac  delta  function. 

We  know  that  there  is  no  boundary  at  x  =  0,  and,  since  the  flow  moves  into 
the  region  of  interest  from  y  =  1,  there  also  can  be  no  boundary  layer  at  y  =  1. 
Thus,  Eq.  (2-12)  satisfies  the  boundary  conditions  n  (0,  y)  »  1  and  n  (x,  1)  =  0  di¬ 


rectly.  At  y  =  0,  we  investigate  the  possible  boundary  structure  by  the  stretching 
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transformation.  * 


y  =  en 


(2-13) 


The  basic  equation 


e  (nxx  +  V  -  nx  +  bny  *  0 


(2-14) 


then  becomes 


(nxx+  72  n™)  *nx  +  7  "n  30 


(2-15) 


or  to  the  lowest  order  in  e 


nnn  +  bnn  *  0 


(2-16) 


The  boundary  conditions  for  this  equation  are,  with  n  as  the  independent  variable 


n  (x,  o)  *  0 


for  x  <  1/b  and  n  (x,n  )  =  0  for  x  >  1/b 


n  (x,  n  >>  0)  =  1 


m 


tx 


& 


A  solution  of  Eq.  (2-16)  that  satisfies  the  boundary  conditions  is 

00 

n  (x,  n  )  3  (1  -  e'bn)  f  -  Six')  dx' 

*  "  b 


(2-17) 


The  two  expressions  (2-12)  and  (2-17)  give  the  density  in  the  bulk  region  and  in  the 
boundary  layer,  respectively.  If,  however,  we  combine  the  two  into 

<k>  k  “ 

nm  J  6  (x')  dx’  -  e  €  f  tf(x')dx'  (2-18) 

y  -  1  +bx  x  -  1/b 

we  have  a  single  expression  which  gives  the  density  approximately  in  the  bulk  re¬ 
gion  and  also  in  the  boundary  layer. 

Although  Eq.  (2-18)  is  an  approximate  solution  on  the  external  boundaries  and 
throughout  the  bulk  domain  of  interest,  there  is  still  the  question  of  formation  of 
internal  boundary  layers,  and  it  is  clear  from  the  nature  of  the  problem  that  such 
an  internal  boundary  layer  will  exist  at  y  *  1  -  bx.  To  investigate  this  free  (in¬ 
ternal)  boundary,  let  s  =  distance  normal  to  y  *  1  -  bx  and  t  =  distance  measured 
from  x  =  0,  y  =  1  along  y  ■  1  -  bx. 

Then  the  transformation  from  (x,  y)  to  (t,  s)  coordinates  is  given  by 


The  variables  may  be  changed  from  x,  y  to  s,  t  in  Eq.  (2-14)  by  using  the 
above  relations  or  by  observing  that  the  line  y  ~  1  -  bx  is  a  characteristic  with 
flow  along  it  and  directly  applying  the  equation  of  continuity.  In  either  case,  one 
obtains 

e  (n  +  n,,)  -  Jl  +  b2  n  -  0  (2-19) 

ss  tt  *  t 


If  we  let  (another  stretching  transformation) 


or 


+  ntt) s  /1  +  fc2  nt 


nnn  +  e2r  ntt  =  e  2r  '  1  /l  +  b2  nt  (2-20) 

Now,  if  in  Eq.  (2-20)  we  set  r  =  1  and  take  terms  of  lowest  order  in  e  ,  we  obtain 
an  equation  the  solution  of  which  does  not  have  acceptable  asymptotic  properties. 
We,  therefore,  let  r  =  1/2  and  obtain  to  the  lowest  order  in  e 

nnn  =  J\  +  b2  nt  (2-21) 

Although  it  is  common  to  obtain  an  ordinary  differential  equation  for  boundary 
layer  solutions,  Eq.  (2-21)  shows  that  this  is  not  invariably  the  case.  The  bound¬ 
ary  conditions  for  Eq.  (2-21)  are 

n  (o,  t)  3  1/2  (fixes  boundary  location) 

lim  n  (n,  t)  =  0 
n  >>o 


lim  n  (n,  t)  =  1 

n  <  <  o 

oo 

n  (n,  o)  3  J  Hr)')  dr,' 
n 

The  last  of  these  is  based  on  the  argument  that  at  t  *  0,  ions  have  just  entered  the 
region  of  interest  and  diffusion  cannot  have  acted  to  establish  a  finite  boundary 
layer.  Noting  that  erf(-z)  =  -erf(z)  and  that  erfc(z)  =  1  -erf(z),  we  find  that  the 
differential  equation  and  all  boundary  conditions  are  satisfied  by 
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V-f >1  ,•  vl 


We  next  observe  that  the  erfc  function  not  only  satisfies  the  free  boundary  layer 
requirements,  but  also  for  arg  (erfc)  <<  0  closely  approximates  the  previous 
bulk  solution.  Further,  the  boundary  layer  equation  at  y  »  0  is  an  ordinary  equa¬ 
tion  in  y  alone:  It  is  still  satisfied  when  multiplied  by  any  arbitrary  function  of  x, 
We  previously  chose  a  delta  function  to  match  to  the  bulk  solution  at  e*  ^  >  1, 

Now  that  we  have  a  solution  that  incorporates  the  free  boundary  properties,  we 
may  just  as  easily  match  the  y*  0  boundary  layer  to  this  bulk  solution.  This  re¬ 
sults  in 

f  bx  -  (1  -  y)  1  --by  C-feS—  j— 

n  =  1/2  erfc  [  2 j  -  e  1/2  erfc  2,/ev£-rb 

as  a  final  solution.  This  solution  is  valid  on  both  boundary  layers  and  throughout 
the  remainder  of  the  domain  of  interest. 

We  wish  to  compare  this  solution  with  the  series  solution,  for  some  larger 
values  of  o  ,  -p.  Since  it  is  fluxes  we  are  mostly  concerned  with,  we  take 


(2-2J 


dn 


dy 


y  =  0 


to  obtain 


-pF 


4) 


erfc 


bx  -  1 

1 

-(bx  -  l)2 

yx  +  t) 

-  4^  exP 

.  4 e (x  +  b>_ 

x  ( 1 


bf  .  _b 

:T  >  ,2. 


(x  +  b) 


3/2 
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For  k*  ~  -  <  <  0,  we  have  -pF  a-  7  =  -p. 
v/c  E 


For  b  =  1,  x  =  1,  and  e  =  l/o  *  1/60, 

-pF  *  28.6;  the  series  value  is  27. 

Figure  2-7  shows  a  comparison  between  fluxes  obtained  from  the  series  solu 
tion  and  from  the  singular  perturbation  solution  for  0  =  -p  =  60.  The  agreement 
is  fair,  and  it  is  expected  that  the  perturbation  values  will  increase  in  accuracy 
with  increasing  0  ,  -p. 

Perturbation  solutions  for  n,  for  o  =  -p  =  100  and  y  =  0.5  have  already  been 
shown  in  Figure  2-6. 


mm  siov  generation.  and  recombination  with 

SPHERICAL  GEOMETRY 

Although  we  have  concluded  in  Section  1  that,  in  general,  diffusion,  flow,  and 
mobility  are  the  dominant  processes  in  stratospheric  sampling,  we  here  solve  a 
problem  that  ignores  the  latter  two  and  includes  two  other  processes  of  lesser  im¬ 
portance.  The  reason  for  this  is  that,  whereas  flow  and/or  mobility  effects  may 
be  completely  absent  in  principle  and  essentially  so  on  occasion,  in  practice,  dif¬ 
fusion  of  charged  particles  cannot  be  absent  when  a  sampling  surface  is  introduced 
into  a  plasma.  Thus,  in  a  sense,  it  underlies  all  other  processes.  The  solution 
of  a  problem  of  this  type  is  thus  informative  and,  as  will  be  seen,  is  not  too  diffi¬ 
cult. 

We  consider  the  following  problem.  A  sphere  collects  charged  particles  by 
diffusion  from  a  plasma  in  which  both  generation  and  square  law  recombination 
occur.  It  is  assumed  that  there  is  no  charge  separation  (hence,  no  electric  fields) 
and  a  single  constant  diffusion  coefficient.  Since  generation  is  by  cosmic  rays  and 
material  bodies  are  for  the  most  part  transparent  to  cosmic  rays,  we  ignore  the 
effect  of  the  sphere  on  generation  and  take  the  generation  rate  to  be  constant  as 
well. 

The  number  density  equation  for  ions  undergoing  production,  recombination, 
and  diffusion  in  the  presence  of  a  spherical  collecting  surface  is 


Here  P  =  production  rate  (number /cm  /sec) 

3 

at  =  recombination  coefficient  (cm  / sec 
r  =  radial  coordinate 
and  D  and  n'  are  as  previously  defined. 

Let  n'  =  n  y  =  \A*//J  y 
x  =  a/r 

...  M  ,  JsiL  a2 

n0D  “  D 

where  a  is  the  sphere  radius. 

Then 


The  use  of  o  with  a  meaning  other  than  that  in  the  previous  section  will  not  lead 
to  confusion. 

The  effect  of  surface  recombination  is  to  fix  the  number  density  of  charged 
particles  at  zero  on  the  surface  of  the  sphere.  The  boundary  conditions  are  thus 

y  (0)  :  i 

y  (1)  =  0  (3-3) 

Before  proceeding  to  solve  Eq.  (3-2),  it  is  worth  knowing,  for  actual  param¬ 
eters,  what  values  of  o  are  likely  to  be  encountered.  If,  for  example,  it  should 
turn  out  that  o  is  less  than,  say,  0.01,  there  is  then  no  need  to  solve  the  full  equa¬ 
tion  for,  in  that  event,  we  have  directly  y  w  1  -  x.  Taking  typical  values  from 
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Appendix  O,  we  find  that  a/a  *  0.014/cm  at  20km  and  4  x  10  /cm  at  40km. 
Since  we  envision  in  solving  this  equation  the  behavior  of  the  whole  sampling  plat¬ 
form,  we  set  a  as  100  cm  to  obtain  a  »  140  at  20  km  and  o  -  0.  4  at  40  km.  It  is 
thus  clear  that  it  is  necessary  to  solve  the  full  Eq.  (3-2). 

Eq.  (3-2)  appears  to  have  a  simple  form,  but  the  non-linear  term  precludes 

solutions  in  terms  of  elementary  functions,  and  the  equation  is  not  listed  in  either 
3  4 

Kamkc  or  Murphy.  We  solve  it  numerically  by  non-linear  relaxation.  Relaxa¬ 
tion  techniques  are  quite  common  in  solving  boundary  value  problems.  However, 
there  do  not  seem  to  be  many  extensive  discussions  of  convergence  of  non-linear 
equations.  Formally,  the  extension  to  this  class  of  problems  is  straightforward. 

In  the  general  case,  solution  of  one  implicit  difference  equation  is  required  for 
each  mesh  point.  For  Eq.  (3-2),  however,  solution  of  implicit  equations  is  not  re¬ 
quired  since  the  difference  equations  can  be  solved  explicitly  by  the  quadratic  for¬ 
mula.  Since  the  solution  is  merely  formal,  however,  even  though  the  results  do 
converge,  there  is  no  guarantee  that: 

(1)  any  solution  to  the  difference  equation  does,  in  fact,  approach  the  solution 
to  the  original  differential  equation; 

(2)  any  solution  found  is  unique. 

We  do  not  explore  either  of  these  two  questions  here.  We  do,  however,  make  a 
minimal  effort  to  insure  that  the  solutions  are  unique  by  a  variation  of  the  initial 
values. 


3.  Kamke,  E.  (1959)  Differentialgleichungen  Losungsmethoden  und  Losungen, 

Chelsea  Publishing  Co.,  New  York. 

4.  Murphy,  George  M.  (1960)  Ordinary  Differential  Equations  and  Their  Solutions, 

D.  Van  Nostrand  Co. ,  Inc.,  New  York. 


The  results  of  these  calculations  are  shown  in  Figure  3-1.  Details  of  the  so¬ 
lution  are  given  in  Appendix  B. 
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Figure  3-1.  Non-Dimensionalized  Number  Density  y  vs  Non-Dimensionalized 
Distance  x.  With  a  as  Parameter  for  a  Spherical  Geometry  Without  Flow 

The  form  of  these  solutions,  as  well  as  of  the  Eq.  (3-2),  suggest  that  solu¬ 
tions  may  be  approximated  by 

y  =  1  -  x  exp[/cT  all  -  l/x>]  (3-4) 

where  a  is  a  constant  to  be  determined.  An  expression  of  this  form  automatically 
satisfies  both  boundary  conditions  and  gives  an  approximate  scaling  with  o .  If 
Eq.  (3-4)  is  put  into  Eq.  (3-2),  there  results 

a2  -  2  +  x  exp[v/aa  (1  -  1/x)]  =  0  (3-5) 

We  thus  have  1  <  a  <  \/2.  We  take  a  to  be  the  geometric  mean  of  these  limits, 
1/4 

a  =  (2)  and  find  that  Eq.  (3-4)  agrees  with  the  numerical  solutions  to  within  3 
percent  for  values  of  0  <  o  <  1000.  The  comparison  is  illustrated  in  Figure  3-2. 


<r 


Figure  3-2.  Comparison  of  Numerically  Calculated  and  Empirical  Formula 
Values  for  Ion  Densities  as  Functions  of  o .  6  =  maximum  difference:  4  - 

maximum  relative  percent  difference;  x^j.  =  value  of  x  corresponding  to  4; 
xM  -  value  of  x  corresponding  to  4  * 


If  the  numerical  values  are  y^  and  the  values  given  by  Eq.  (3-4)  y^,.  then  the 
graphs  show 

6  =  ^C  ’  yF)  MAX 

and 

)*100 
MAX 

as  functions  of  o.  Also  shown  are  the  corresponding  values  of  x,  x^  .  and  xM^. 

Of  prime  interest  as  far  as  collection  of  ions  is  concerned  is  the  flux  to  the 
surface.  Because  of  the  accuracy  of  the  approximation  in  Eq.  (3-4),  we  risk  the 
dangers  associated  with  differentiating  approximations  and  use  this  formula  to  ob¬ 
tain  the  flux.  For  this  purpose,  we  convert  all  quantities  except  o  back  to  dimen¬ 
sional  variables  to  obtain  for  the  flux 
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The  flux  to  the  surface  is  thus  increased  by  the  ratio  (1  +  </a )  over  what  it  would 

be  in  the  absence  of  production  and  recombination  near  the  surface. 
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Appendix  A 

Numerical  Solution  of  Diffusion,  Mobility,  Flow  Equation 

for  Rectangular  Geometry 


This  is  an  elliptic  equation  which  we  will  solve  by  a  relaxation  (iterative)  method. 

Numerical  values  of  a  are  chosen  so  that  a  *  O  [(P/o  )].  This  strikes  a  balance 

between  two  extremes  related  to  the  value  of  n  *  n.  .  at  the  first  grid  point  adja- 

cent  to  z  *  0.  In  one  extreme,  nj  ..  is  so  small  as  to  be  of  no  interest;  in  the  other, 

n.  .  is  so  large  that  information  is  lost, 

*»  J 

We  next  set  up  a  grid  as  shown  in  Figure  A  - 1  which  also  shows  values  of  n  on 
the  boundaries. 

The  most  common  iteration  scheme  for  solution  of  partial  differential  equa¬ 
tions  uses  centered  differences  for  all  derivatives.  For  Eq.  (A -2),  this  results 
in 
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>  > 


X  =  ®  x  =  0 

i  *  0  z  =  ih  =  i/k  i  =  k 


Figure  A-l.  Calculational  Grid  Mesh  for  Numerical  Solution  of  Flow/Diffusion 
Equation  for  a  Rectangular  Geometry;  Mesh  Has  20-by-10  Increments 


n.  .  =  1/2  [i2a2  +  l2 


i.j 


1  ' 


«■*»>)  *"i.  j»i(“2‘“)*  V  j-i(«2*-f )l  IA-31 


with 


x.  =  -(1/a)  In  Zj  -  -  1/a  In  (i/k) 


This  formulation  is  suspect  on  two  grounds.  First,  noting  that  p<  0,  the  two 

2 

(o  +  o  k  ana  n.  .  .  <  «- 

2 

large  o  and  -p.  be  of  the  wrong  sign  to  guarantee  the  existence  of  unique  solu- 


,.2  2 

terms  n.  .  .  (l  a  -  —  (a  +  a  ))  and  n.  .  ,  (l 

i-l  -J  2  ‘'I'1 


pi! 2)  will,  for  sufficiently 


tions  (Ames,  p.  101).  Next,  suppose  we  calculate  on  the  first  pass  the  point 


ni,  1  *  ni.  j-  1 


i  *  k  -  1,  j  -  1.  Then  from  the  boundary  conditions,  n.  .  . 

t  1 » .1 

=  0  and  n.  ,  .  =  1.  From  this,  n.  . **o  and  n.  ,  can  have  anv  value,  no  matter 

i +  1.]  i.  3  i.J 

how  large,  depending  on  how  large  we  make  o.  Although  possibly  not  fatal,  this 


does  not  seem  like  a  good  prognosis  for  stability  and  convergence. 


! 

L 


C. 


% 


If,  on  the  other  hand,  we  use  centered  second  differences  but  forward  first 
differences, 

_£U_  _♦  ih  ni+l,j  ~  ni,j  :  gn  —  n*.  j-H  '  ni.J 
dz  h  dy  P 

we  obtain 

n.  .  *  (2  i2a2  +  2*2  +  ia  (a  +  o  )  -  pi  )  * 

1‘3 

jni+l,j<i2a2  +  ia(“+a>  +  ni-l.ji2“2 

+  ni.j+l(*2-PM  +  ni, j-i  *2 1  (A'4) 

The  off  diagonal  terms  now  all  have  the  correct  sign  and  nj  ^  remains  of  the  order 
unity  no  matter  how  large  o  and  -p  become.  We  used  Eq.  (A -4)  for  all  numerical 
calculations. 


Appendix  B 

Numerical  Solution  of  Diffusion,  Generation.  Recombination 

Equation  for  Spherical  Geometry 


Eq.  (3-2)  may  be  represented  by  the  following  (centered)  difference  equation: 


J  V  ( Ax >2  /  i 


Setting  and  solving  for  y^ 


1) 


yj  =  |  *'*4  (4x)2  *  J8  <4x>4  +  o  (°  +  j4  <4x)2  (y^j  +  |  <B'2) 


Eq.  (B-2)  holds  for  y .  i  y  ,  y  ,  the  two  boundary  values.  For  these  values,  we 
j  o  m 

have  y  =  1  and  y_  *  0.  The*  sign  in  Eq.  (B-2)  is  chosen  +  to  make  y,  >  0. 
o  m  j 

To  give  some  indication  of  convergence  and  uniqueness,  we  solve  for  each 
value  of  o  twice  (superscripts  are  iteration  numbers),  once  with 


y°j s  j  * x  0 

and  again  with 

■  »•  j  *  »•  ■  >• 
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In  all  cases,  the  solutions  for  these  two  initial  conditions  continued  to  approach 
each  other  more  closely  with  successive  iterations. 

Typical  values  for  m  were  m  »  25  and  convergence  criteria 


MAXj  were  typically  ~10 


Typically,  on  the  order  of  several  hundred  iterations  were  needed  to  attain 
this  convergence. 


Appendix  C 

Typical  Values  for  Parameters  Occurring  in  Text 


Parameter 

Description 

20  km 

40  km 

n  (1/cm3) 

o 

Equilibrium  ion  density 

5500 

4000 

3 

fl(cm  / sec) 

Recombination  coefficient 

5  x  10"7 

5  x  10~8 

/S(l/cm3sec) 

Ion  Production  Rate 

15 

0.  75 

3 

l>(cm  /sec) 

Kinematic  Viscosity  m/p) 

2 

50 

Vg  (cm/sec) 

Platform  Drift  Velocity 

100 

100 

Re 

Drift  Reynolds  Number 

50L 

2L 

6  (cm) 

Drift  Boundary  Layer  Thickness 

JZn.o 

/L/1.4 

Po'P 

Sea  Level  to  Ambient  Density  Ratio 

14 

300 

/.(cm) 

Neutral-neutral  Mean  Free  Path 

9x  10'5 

2x  10'3 

L  (cm) 


Characteristic  Dimension 


100  for  platform 


